%% Parameters - main parameters to calibrate the model
d_hat = 40; % Radius of influence for a 100 mmgy plant
theta = 0.5; % Exponent on capacity used to calculate d_bar = alpha*K^theta

%% Load data
% Load Excel data with ethanol plant capacities and locations
ethanolPlants = readtable("Ethanol Plant Data/Ethanol Plant_location cap capCF.csv");

% Extract important variables into matrices
% Plant capacities by year
years = 2001:2018;
vars= "cap"+years;
varsCF = "capCF"+years;
cap1 = table2array(ethanolPlants(:,vars));
cap1(isnan(cap1)) = 0;
capCF1 = table2array(ethanolPlants(:,varsCF));
capCF1(isnan(capCF1)) = 0;
% Plant location - [lat, long]
locations = [ethanolPlants.lat ethanolPlants.lng];

%% Calculate radius of influence d_bar for each plant in each year
% d_bar = alpha*K^theta
alpha = d_hat/(100^theta); % Calibrate alpha based on d_hat
d_bar1 = alpha*cap1.^theta;
d_barCF1 = alpha*capCF1.^theta;

%% Calculate distance matrix between all plants
% Radius of the earth in miles; needed for distance calc
R = 3959;

% Define constants for the haversine calculation
phi1 = locations(:,1)*pi()/180; % Latitude in radians
phi2 = phi1'; % Latitude in radians, change to row vector
deltaPhi = phi1-phi2; % Latitude difference matrix
deltaLambda = (locations(:,2)-locations(:,2)')*pi()/180; % Longitude difference matrix

% Haversine calculation
a = sin(deltaPhi/2).*sin(deltaPhi/2) + cos(phi1).*cos(phi2).*sin(deltaLambda/2).*sin(deltaLambda/2);
c = 2*atan2(sqrt(a),sqrt(1-a));

distances = R*c; % Distance measured in miles

%% Add up all effective capacities within d_bar1 of each plant to get d_bar2
ph = zeros(height(ethanolPlants),length(years)); % Placeholder matrix of zeros
near_plants_sum = ph;
cap2 = ph;
capCF2 = ph;

for y=1:length(years)
    % Compute acutal effective capacities
    d_bar_sum = d_bar1(:,y) + d_bar1(:,y)';
    near_plants = (distances<d_bar_sum).*(d_bar1(:,y)>0);
    near_plants_sum(:,y) = sum(near_plants,2);

    cap2(:,y) = (near_plants.*(1-distances./d_bar_sum))*cap1(:,y);

    % Compute counterfactual effective capacities
    d_barCF_sum = d_barCF1(:,y) + d_barCF1(:,y)';
    near_plantsCF = (distances<d_barCF_sum).*(d_barCF1(:,y)>0);

    capCF2(:,y) = (near_plantsCF.*(1-distances./d_barCF_sum))*capCF1(:,y);
end

cap2(isnan(cap2)) = 0;
d_bar2 = alpha*cap2.^theta;

capCF2(isnan(capCF2)) = 0;
d_barCF2 = alpha*capCF2.^theta;

%% Export data
% Export:
%   d_bar1: radius without considering interference from other plants -
%       OVER estimate of true radius
%   d_bar1CF: same as above, but for counterfactual ethanol timeline
%   d_bar2: radius accounting for interference from other plants - UNDER
%       estimate of true radius
%   d_bar2CF: same as above, but for counterfactual ethanol timeline
           
varNames1 = "d1_"+years;
varNames2 = "dCF1_"+years;
varNames3 = "d2_"+years;
varNames4 = "dCF2_"+years;

T1 = array2table(d_bar1,"VariableNames",varNames1);
T2 = array2table(d_barCF1,"VariableNames",varNames2);
T3 = array2table(d_bar2,"VariableNames",varNames3);
T4 = array2table(d_barCF2,"VariableNames",varNames4);

exportTable = [ethanolPlants T1 T2 T3 T4];
writetable(exportTable,"Plant Effective Radius.xlsx");